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Abstract The research about asteroids attracts more and more attention recently, es- 
pecially focusing on their physical structures, such as the spin axis, the rotation period 
and the shape. The long distance between Earth observers and asteroids makes it im- 
possible to get the shape and other parameters of asteroids directly with the exception 
of the NEAs (Near Earth Asteroids) and others passed by some spacecrafts. Generally 
photometric measurement is still the main way to obtain the research data for asteroids 
now, i.e. the lightcurves recording the brightness and positions of asteroids. Supposing 
that the shape of the asteroid is a triaxial ellipsoid with a stable spinning status, a new 
method is present in this article to reconstruct the shape models of asteroids from the 
lightcurves, with the other physical parameters together. By applying a special curvature 
function, the method calculates the brightness integration on a unit sphere and Lebedev 
Quadrature is employed for the discretization. At last the method searches the optimal so- 
lution by Levenberg-Marquardt algorithm to minimize the residual of the brightness. By 
adopting this method not only related physical parameters of asteroids can be obtained at 
a reasonable accuracy, but also a simple shape model of Ellipsoid can be generated for 
reconstructing more sophisticated shape model further. 
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1 INTRODUCTION 

With the development of the technology humans have much more ability to discover the universe now. 
The beginning of our solar system is still an interesting problem, which results in the research about 
asteroids because the asteroids reserve much information about the original formation of the solar sys- 
tem. Most of asteroids found now lie between Mars and Jupiter, commonly called as main belt aster- 
oids(MBAs). Generally MBAs are orbiting the sun with a stable status with the exception that some 
larger bodies pass them over and the perturbations will change their orbits. Due to the long distance of 
2.1 ~ 3.3AU between the earth and most of asteroids, until now the main research data about asteroids 
is the photometric brightness with the position coordinates of asteroids, the earth and the sun, generally 
i.e. lightcurves, which are recorded by many ground-based observatories. 

Shape models with other physical parameters such as the rotation period and spin axis, especially 
the albedo of asteroids are the most concerned about asteroids. The total distribution of the spin axis and 
rotation period of MBAs can help us research the origin and evolution of the solar system. Commonly 
the rotation periods of MBAs are stable with some hours less than one day. But recently it is found that 
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the thermal r adiati o n from the sun can change the asteroid's rate of rotation, called as YORP effect. 
lLowrv et al.l d2007l) : pTavlor et ail d2007l) observed the rotation period of asteroid 54509 will double in 
about 600,000 years which confirms that the periods of asteroids may change in a very small amplitude 
within a long term. 

There are many methods to reconstruct the shape models of asteroids with the other physical param- 
eters. [RussejJ dHOQ) firstly started the research about the shape of asteroids and concluded a pessimistic 
result that the shape could not be obtained merely depending on the lightcurves of the opposition. With 
the development of both mathematical theories and technologies in telescope more and more lightcurves 
viewed in various solar phase angles are recorded by many observatories located everywhere on the 
earth and the shape inv ersion algorithms can figure out the shape model in a reasonable accuracy by the 
high-speed computers. iLumme & Bowelll (Il981albl) presented an estimation about the scattering law in 
the surface of atmosphereless bodies and introdu ced a spherical harmonics method for asteroid pole 
determination (ILumme et al.L 1 1 9901) . iHapka d 1 9841) moved further about the scattering law with the con- 
sideration of bidirectio nal reflectance an d the roughness of the surfa ce. Basing on the scattering law 
of Lumme and Bowell, iKaittunenl (fl989h : iKarttunen & Bowelll Jl989t) presented a method to generate 
lightcurves of triaxial ellipsoid models and discussed the variations of convex model and nonconvex 
model. But they just noted that it was possible to determine the axial rati os of a triaxia l ellips oid model, 
while did not give a method to find the shape models from lightcurves. ICellino et al.1 (1 1989b adopted a 
model formed by merging together eight octants of ellipsoids having different semi-axes and generated 
the corresponding lightcurves wit hout an inverse meth od, too. Adopting the sparse photometric data col- 
lected by the Hipparcos satellite. [Cellino et al.l (120091) presented a genetic inversion to find the physical 
properties of asteroids with the assumption th at the shape model of the asteroid is a triaxial ellipsoid. 
Further more lKaasalainen & Lambergl(ll992allbl) built up a very efficient method to reconstruct the arbi- 
trary surface of asteroids and the inverse shape models are confirmed by the flyby observation in space 
dKaasalainen et al.ll2002T) . 

Nevertheless until now the ellipsoid model plays an important role. Firstly the photometric data 
observed by the observatories based on the earth is not enough accurate because of the atmosphere mist 
, the CCD heat, wrong operations and so on. Generally the error of photometric data is about 2% and we 
can not expect to obtain a very accurate shape model only from the lightcurves. Secondly the ellipsoid 
model is simple but it can make the inversion easy and efficient with an acceptable physical parameters 
while the statistical research about the spin axis and period of asteroids can obtain the needed data 
in a fast way. Thirdly the rough model inferred from the ellipsoid model may be the initial value for 
reconstructing more accurate shape model by Kaasalainen's method. 

Under the similar definition employed in Kaasalainen's method and tiling the triangular facets in 
Lebedev way we present a fast method of the ellipsoid model in this article as the following organization. 
In section|2]we describe the technical details of the fast ellipsoid model, including the scattering law in 
12.11 and the photometric integration by applying the curvature function in 12.21 Furthermore the inverse 
problem with a discretized integration is presented in 12.31 Finally the formulas of derivatives are given 
out for reference in optimization part 12.41 At last the summing up and future plan will be discussed in 
section |3] 

2 ELLIPSOID MODEL FOR ASTEROIDS 

The disk-integrated photometric data i.e. brightness in lightcurves, contains much information about the 
asteroid which can be applied to obtain the related parameters. The periodic variation of brightness is 
mainly due to the variation of the shape as the asteroid spins around its axis. We assume that the shape of 
asteroid is a triaxial ellipsoid with three semi-axis a > b > c > 0, which spins in the period P expressed 
in hours around its shortest axis whose spherical coordinate is denoted as (A,f}) in the J2000 ecliptic 
frame system. As the convention the brightness data recorded in lightcurves is reduced to unit distances 
of between the aster oid and the sun, the earth respectively and corrected according to the light-time 
(IDurech et al So we suppose lightcurves is processed in this way before applying our method. 
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2.1 The Scattering Law 

The scattering behavior is an inevitable problem in the asteroid model. iLumme & Bowelll d!981allbh 
considered several physical parameters, such as the single-scattering albedo D.q, the asymmetry factor g, 
the volume density of the surface material D, the roughness of the surface g etc. At last a sophisti cated 
scattering model was built to simulate the reflection behavior of the light from Sun. iHaokd (11984 took 
into account the opposition effect and the shadow in the particles of the surface. These scattering models 
can express the physical characteristics of the light refle ction in a rational mann er, but they are not 
efficient in reality due to the uncertain physical parameters. iKaasalainen et all (120051) made a photometry 
research about an artificial asteroid in the laboratory experiments and confirm that the shape variation is 
the main cause for the variation of brightness, not the scattering law. Further more, they found that it is 
hard to distinguish the difference between the scattering law and the random error. In order to acquire the 
shape model in an efficient way, a simple scattering law is needed. IKaasalainen et al.l (1200 lb presented 
an convenient method to simulate the scattering behavior by merging both the single scattering factor 
5Ls(Lommel-Seeliger) and the multiple scattering factor S /.(Lambert) in a linear combination. The 
scattering law can be expressed in 

S(n,no,a) = f(a)[SLs(M,fio) + cS L (jx,(Xo)] 

= /(4^+ rW ), (1) 

where the fi, /vo are defined as follows under the definition of 77(7?, if) as the outward unit normal vector of 
the surface and u>, cog as the directions to the earth and the sun observed from the asteroid respectively, 

11 = 10-1], fi = a> ■ i]. (2) 

The phase function f(a) is a fitted function in the three-parameter form 

/( a )=A exp(— ) + (3) 

where Aq and D are the amplitude and scale length of the opposition effect and k is the overall slope 
of the phase curve. The above scattering law with four parameters adopted in this article can perform 
efficiently in the shape inverse problem and simulate the opposition effect rationally. 



2.2 The Photometric Integration 

In order to reconstruct the shape model of asteroids, the synthetic brightness must be described in the 
direct problem. As mentioned above, the scattering law with four parameters can be adopted herein to 
generate the photometric brightness as a surface integral 

L(cu Q , to) = S (//, no, a) da, (4) 

where E+ is the part of the asteroid shape which is both illuminated by the sun and visible from the earth, 
i.e fi,/jQ > 0. As assuming that the shape model of asteroids is triaxial ellipsoid, the integral (0]i can be 
calculated numerically by the traditional triangulation, tiling the approximately equal triangular facets 
on the surface, which is a linear algorithm as the number of tessellated facets. In Fig [2] it is shown that 
an error level of 10~ 4 needs m ore than 10 4 triangular facets tessellated on the surface of the ellipsoid. 

iLebedev & Laikovld 19991) presented a fa st method to calcul a te the surface integral on the unit sphere 
S by tiling the triangular facets not equally. IKaasalainen et al.1 (|2012|) applied this technique into their 
method and confirmed the lebedev quadrature is efficient in the surface integration. The different distri- 
butions of triangular facets of the traditional triangularization and lebedev quadrature are shown in Fig 

m 
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Fig. 1 The comparison of two different tessellation methods (Left: Triangulation Right: 
Lebedev) 



Due to that the lebedev quadrature is based on the unit sphere, in our method a curvature function 
from the surface of ellipsoid E to the unit sphere S is built with the format 

J dr ^ dr J 

(Ke,4>)= . ^ ,r(6,<f>)eE. (5) 
suit) 

With the curvature function G{9, <p) the brightness integral © can be transformed to the surface integral 
on the unit sphere S 

L(a>o,a>) = I I S(fi,iJ.o,a)G(9,(/>)dcr, (6) 



ff, 



where 5+ is the similar part of unit sphere with fi,/Jo > 0. 

The curvature function under the standard parametrization of the ellipsoidal surface E 

x\ - a sin 9 cos <f>, xz = b sin 9 sin <p, X3 = ccos.9,9 € [O,27r],0 e [0,7r], (7) 

has the form 

I sin 9 cos d> \ /sin0sin< 



/sinfc'cosffl\ /sin&sin(*\ /cosfc'\" 
Gi^-atc^-—^.) +(— 5-^) +(— ). 

while iKaasalainen & LamberS dl992bb also presented the other curvature function with the form 



abc 



^ \ (a sin & cos (p) 2 + (b sin 1? sin ip) 1 + (c cos ff) 2 ) 

under the definition of (1?, ip) denoting the spherical coordinates of the normal vectors of the ellipsoidal 
surface. 

Both of two curvature functions can perform well in computing the brightness integral ©. But the 
first one mentioned and adopted in our method is more efficient for the extreme situations, such as the 
ellipsoid model with a large difference between its three semi axes. The 'Do g-Bone' shaped asteroi d 
(216) Kleopatra has been found that its shape is close to an elongated ellipsoid dDescamps et all 1201 ll) . 
Fig|2]shows the comparison of our curvature function(V) and Kaasalainen's curvature function('o') to 
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compute the surface area of ellipsoids with different semi axes with the following formula easily letting 
the scattering function S (jj.,(1q, a) = 1 in the integral|6] 

Surface Area of Ellipsoids — JJ" G(8,(f>) dcr. (10) 

Apparently our curvature can obtain lower error lever than Kaasalainen's curvature, especially in 
the case of elongated ellipsoids. Besides, the difference between the triangulation('O') and lebedev 
quadrature(V or 'o') is also compared in Fig [2] The dominant lebedev quadrature can guarantee the 
efficiency of our fast method. 
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Fig. 2 The comparison of Triangulation('O') and Lebedev Quadratures with our 
Curvature(V) and Kaasalainen's Curvature('o') respectively for computing the surface area 
of ellipsoids with various semi axes (Left: a=8,b=7,c=6 Right: a=10,b=2,c=1.5) 



2.3 The Inverse Problem 

With the definitions of both the scattering function and the brightness integration the inverse problem of 
the shape model of asteroids can be illustrated totally as follows. 

rrt 1 

Supposing that all lightcurves in this article are formatted in the same way of DAMIT ( Du rech et all 
I2010I) . the relative brightness (L) with respect to a Julian Day (r) and the positions of the earth and 
sun (u>,u>o) in ecliptic asteroid-centric coordinate system can be obtained directly from lightcurves. 
According to the aforementioned brightness integration, all required parameters are three semi axes 
(a,b,c), the spherical coordinate (A, ft) of the spin axis in ecliptic system, the period (P), the four pa- 
rameters (Ao, D, k,y) of the scattering function and the initial phase angle (<5o) at the beginning epoch 
(to). So there are 1 1 parameters in all as fo is generally set to be the first Julian date of all lightcurves. 

The observed brightness will vary periodically as the rotation of the asteroid around its spin axis, 
which seems that the sun and earth spin around the stationary asteroid. Assuming the asteroid coordi- 
nate system is the Cartesian frame with the spin axis as z-axis and the long semi-axis as x-axis, the 
origins of ecliptic and asteroid coordinate system are located at the same point of the center of the 
asteroid. Denoting a>, (Sq as the coordinate of the earth and sun in the asteroid coordinate system, the 
transformation between ecliptic system and asteroid system can be given as 



to = R 3 (®)R 2 (J3)R 3 (A)a), (3 = R 3 (®)R 2 (J3)R3(A)a>o, 



(11) 
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where Ro, R3 are the rotation matrices rotated as the y and z axis respectively with the forms 



Ri{0) = 



1 
sin# cosfl , 



cos 9 sin 9 ' 
- sin cos 9 
1, 



(12) 



As the asteroid rotates in an angular speed 2tt/P, phase angle <t> varies with respect to the epoch t (Julian 
Day) in the form 



In 

3>(f) = $0 + y{t - to), 
while if the YORP effect is taken into account, the corresponding form will be 

$(0 = $0 + y (t - to) + -v(t - t )\ 



(13) 



(14) 



where v = dtl/dt (Q = 2tt/P) is the rotation rate of change of the period P. 

Let (9j, <pi) be the discretized points of Lebedev quadrature with the total number of facets N, tes- 
sellating on the surface of unit sphere with a corresponding weight w,, which can be also treated as the 
area of small facets. With the normal vector of any small facet of the ellipsoid shape 



m = 



sin 9j cos d>i sin 9-, sin <pj cos 9-, 



a 



(15) 



the brightness integral © can be discretized as 

N 

L(co, ojo) « J] (S 0/°,//®, a)G(9 h <pdw) , 



(16) 



!=1 



where are the inner products of u>, a»o and the unit normal vector i/,7|//,i respectively. Further 

more under the definition of normal vector (TT3T >. the curvature function (J8j can be simplified in the form 



G(6i,<p t ) = abc\ijii 



(17) 



Finally merging the scattering function S (/i,/io, a) in the formula (HJ the discretized brightness integral 
(|T6) will be 



L(w, w ) w ^ /(«)(&» • i}i)(^o ■ »7i) 

i=i \ 



1 



y 



(a) ■ 77,-) + (wo • go \m\ 



abcwi , 



(18) 



where (jc • y) is the inner product of two vectors x,y. 

The inverse problem of the ellipsoid model can be described to find the 1 1 parameters mentioned 
above to minimize the 



X 



-s 



I.® 




<L©> 


(L©) 



(19) 



where L ( '\ L (,) denote the brightness vectors of observed data and the synthetic data containing all points 
in the zth lightcurve, while (L (,) ), (L w ) denote the mean brightness of the two brightness vectors. 
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2.4 The Optimization Algorithm 

There are many methods to find the best fit solution of the inverse problem ( fT9l ), such as the genetic algo- 
rithm, methods of conjugate gradients and so on. Herein we employ the Levenberg-Marquardt method 
to search the optimal solution, whi ch works very we ll in practice and has almost become the standard 
of nonlinear least-squares routines dMarquardtl 1 1 963b . The LM(Levenberg-Marquardt) method can con- 
verge in a very fast speed, while the obtained result is often the local optimal solution. There is a general 
method to complement its deficiency, i.e. searching the best fit result with various initial values. The 
details of LM method can be found in iPress et al.ld 1992b . 

But it needs to compute the derivatives of x 1 m LM method, which means the derivatives of the 
discretized brightness integral ( fTST l with respect to the 1 1 parameters have to be computed in formula! 
In order to make it program easily we give out the formulas for reference. 

Letting 



fi + flQ 



yjm 
\m\ " 



(20) 



the derivative of L(a>, coo) with respect to long semi axis(a) has the form 



AT N 

~d~a 

/=1 



dS dyu dS d/i y^o sin 2 6» ; cos 2 . 

f(a)bcWiS + abcwi \^-^~ + t ^— + t— ^ 

dfi aa djio oa a |//;| 



(21) 



The other semi axes (such as b, c) have the similar formulas. The derivatives of L(co, coo) with respect to 
(A,p,P,<S>o) are 



3A ~ 



dS dR C/l") dS dR ( '/0 \ 

>i[—R3(®)R2(P)—^-CO ■ m + —R 3 (^)R 2 (fi)—^-loj ■ T]ij 



(22) 



8L 



(23) 



_ - V 



f(a)abcwj 



dS dR 3 (<&) 



R 2 (/3)R3(A)to ■ m 



dS 8R 3 m 



s 3// SO ^ ' d/ao SO 

Due to the relation of P and (Do in one function (1131) . 

dL dL d<S> 



R 2 (J3)Ri(A)co Q ■ rji 



dP 5(D dP 



(24) 



(25) 



The derivatives of L(co, coq) with respect to the three parameters Aq, D, k in scattering function are easy 
to obtain by multiplying the corresponding derivatives of f{a) respectively, 



df a df A a a df 



(26) 



and the derivative to y is 



dL 
dy 



f(a)abc\Vj 



\m\ 



(27) 
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3 CONCLUSIONS 

We have described a fast method to obtain the physical parameters and the shape models of asteroids 
basing on the ellipsoidal shape. This method adopted the Lebedev quadrature to discretize the surface 
integral on the unit sphere, which can decrease the computational cost largely with a high accuracy. In 
addition now that this method can compute the period and the orientation of the spin axis in an efficient 
way, we can rely on its result to refine the shape of asteroids, such as adopting Kaasalainen's method. 

The related formulas are presented in this article totally, such as curvature function (O and dis- 
cretized brightness integral ( TT~8T >. especially the complicated derivatives of^- 2 . Here it should be noted 
that all the formulas are given out in a general form with three semi-axis parameters a, b, c because this 
format is comprehensible and easy to compare with other formula like Kaasalainen's Curvature func- 
tion. In computation due to the relative brightness the semi-axis parameters a, b, c are not exactly same 
as the real asteroids in size. So the ratios of axis of ellipsoid model a/b, a/c are more practical. By ap- 
plying the fast method, it is easy to compute the ratios from the derived parameters at last or we can set 
c = 1 in Levenberg-Marquardt method to obtain the ratios directly. Besides, we can also replace the 22 
in <D(f) by its angular speed Q to simplify the computation. The detailed numerical test and application 
to real asteroids will be presented in the future article. 
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